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Abstract 

In the presented article we present an algorithm for the computation of ground state spin configurations for the 2d random bond 
Ising model on planar triangular lattice graphs. Therefore, it is explained how the respective ground state problem can be mapped 
to an auxiliary minimum-weight perfect matching problem, solvable in polynomial time. Consequently, the ground state properties 
as well as minimum-energy domain wall (MEDW) excitations for very large 2ii systems, e.g. lattice graphs with up to A^ = 384x384 
spins, can be analyzed very fast. 

Here, we investigate the critical behavior of the corresponding T -Q ferromagnet to spin-glass transition, signaled by a breakdown 
,of the magnetization, using finite-size scaling analyses of the magnetization and MEDW excitation energy and we contrast our 
■numerical results with previous simulations and presumably exact results. 

'fCeywords: Random bond Ising model. Negative- weight percolation, Groundstate phase transitions 
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.1. Introduction 

' Triggered by the exchange of ideas between computer sci- 
ence and theoretical physics in the past decades, it was realized 
that several basic problems in the context of disordered systems 
relate to "easy" optimization problems. These are problems 
where the solution time is polynomial in the size of the problem 
description. As a result, many disordered systems can now be 
analyzed numerically exact through computer simulations by 
using fast combinatorial optimization algorithms 111 |2|, |3|] . E.g., 
ground state (GS) spin configurations for the random-field Ising 
magnet (in any dimension d) can be obtained by computing the 
maximum flow for an auxiliary network problem [2]. Another 
example is the 2d Ising spin glass (ISG), where the lattice can 
be embedded in a plane. For this model, the problem of find- 
ing a GS spin configuration for a given realization of the nearest 
neighbor couplings can be mapped to an appropriate minimum- 
weight perfect-matching (MWPM) problem lJ,|3l- Finally, the 
MWPM problem can be solved in polynomial time by means 
of exact combinatorial optimization algorithms |5]. Thus, the 
planar 2d ISG can be studied directly at zero temperature with- 
out equilibration problems and within polynomial time. Hence, 
very large systems can be considered, giving very precise and 
reliable estimates for the observables. Actually, there are dif- 
ferent approaches that allow for an exact computation of GSs 
for the planar 2d ISG [6|,LZI|8|]. Albeit all of these approaches 
rely on the computation of MWPMs on an auxiliary graph, they 
differ regarding the subtleties of the mapping to the respective 
auxiliary problem. The most efficient of these approaches (see 
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Ref. ISl]) is based on the Kasteleyn treatment of the Ising model 
JSt], which previously was also used to obtain extended ground 
states for the 2d ISG with fully periodic boundary conditions 

Here, we introduce a dedicated algorithm that yields exact 
GS spin configurations for the 2d random-bond Ising model 
(RBIM) on planar triangular lattice graphs. As the previous 
approaches, the algorithm presented here requires to solve an 
associated MWPM problem. The corresponding mapping uses 
a relation between perfect matchings and paths on a graph 
I 111 Il2r . In effect, these paths can be used to partition the 
graph into domains of up and down spins that comprise a GS 
spin configuration, see Fig.[T] Consequently, the GS properties 
as well as minimum-energy domain wall (MEDW) excitations, 
see Fig.[T] can be analyzed very fast. The presented algorithm 
enables us to study large systems, while allowing for an appro- 
priate disorder average within a reasonable amount of comput- 
ing time. In this regard, it requires to compute a MWPM for an 
auxiliary graph with 0{N) edges only (wherein A^ is the num- 
ber of spins on the lattice). However, note that the algorithm 
presented here is asymptotically not faster than the algorithm 
presented in Refs. Isl llOll . but it highlights the algorithmic rela- 
tion between the GS problem for spin glasses and the recently 
proposed negative-weight percolation (NWP) problem 11 1311 . 

In the presented article, we investigate the critical behavior 
of the T -Q ferromagnet (FM) to spin-glass (SG) transition for 
the 2d RBIM, signaled by a breakdown of the magnetization, 
using finite-size scaling (ESS) analyses of the MEDW excita- 
tion energy. In this regard, we obtain a highly precise estimate 
of the critical point for the triangular lattice geometry and we 
verify the critical exponents obtained earlier for the RBIM on 
the planar square lattice il4lll5ll . Finally, we contrast our nu- 
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Figure 1 : Samples of a ±J random bond Ising spin system on a triangular 
lattice of side length L = 32. The samples are taken at three different values of 
the disorder parameter p. (a) p = 0.15 characterized by a ferromagnetic GS, (b) 
/7 = 0.2 characterized by a GS with SG order and (c) p = 0.5, i.e. the canonical 
±J ISG. In the figure, periodic BCs are indicated by the dashed vertical lines. 
From left to right: (Left) Transition graph that describes the diiference between 
a ferromagnetic reference spin configuration and the GS spin configuration, 
(center) corresponding GS, (right) MEDW excitation relative to the GS. 



merical results with previous simulations and presumably exact 
results ifiell . 

The remainder of the presented article is organized as fol- 
lows. In section |2l we introduce the model in more detail and 
we outline the algorithm used to compute the GS spin config- 
urations. In section [3] we present the results of our numerical 
simulations and in section|4]we conclude with a summary. 



2. Model and Algorithm 

In the presented article we perform GS calculations for the 2d 
RBIM, where the respective model consists of N - Lx L Ising 
spins cr - (cri, . . . ,crf^), where cti - ±1, located on the sites 
of a planar triangular lattice graph. Therein, the energy of a 
given spin configuration is measured by the Edwards- Anderson 
Hamiltonian 



H{cr) = - ^ Jij o-iO-j, 



(1) 



<'■.;> 



where the sum runs over all pairs of nearest-neighbor spins (on 
the triangular lattice) with periodic boundary conditions (BCs) 
in the x-direction and free BCs in the y-direction. In the above 
energy function, the bonds 7,j are quenched random variables 
drawn from the disorder distribution 



P(J) 



p6(J+l) + (l-p)5{J-l). 



(2) 



Therein, one realization of the disorder consists of a random 
fraction p of antiferromagnetic bonds (J = -1) that prefer an 
antiparallel alignment of the coupled spins, and a fraction (l-p) 



of ferromagnetic bonds (J - 1) in favor of parallel aligned 
spins. In general, the competitive nature of these interactions 
gives rise to frustration. A plaquette, i.e. an elementary trian- 
gle on the lattice, is said to be frustrated if it is bordered by an 
odd number of antiferromagnetic bonds. In effect, frustration 
rules out a GS (i.e. a minimizer ctqs of Eq.[T]i in which all the 
bonds are satisfied. As limiting cases one can identify the Ising 
ferromagnet at ;? = and the canonical +J ISG at p - 1/2. 
Hence, as a function of the disorder parameter p we expect to 
find a ferromagnetic phase (spin-glass phase) for p<pc(p> Pc), 
wherein pc denotes the critical point at which the T -0 FM-SG 
transition takes place. For the ISG with a bimodal disorder dis- 
tribution the GS is highly degenerate and the average number of 
such GSs increases exponentially with A^ BlTlllSll . Apart from 
the GSs, we here also aim to characterize the energetic proper- 
ties of MEDW excitations. A domain wall is an interface that 
spans the system in the direction with thee free BCs. Now, the 
MEDW is such an interface with an excitation energy SE that 
is minimal among all possible domain walls. Due to the exten- 
sive degeneracy of the GSs, the "lattice-path" associated with a 
MEDW is not unique. I.e., there are many DWs with minimal 
excitation energy. Albeit the geometric properties of a MEDW 
are not unique, its excitation energy is unique. MEDWs for 
three different values of the disorder parameter p are illustrated 
in Fig. [1] 

We now give a brief description of the algorithm that we use 
to compute the GSs. Therefore, we first set a reference spin 
configuration ctr. The most convenient choice is a maximally 
polarized, i.e. ferromagnetic, configuration ctr - (+1, . . . , -i-l). 
Then, we construct a weighted dual of the spin lattice as shown 
in Fig. |2|a). Since the spin lattice considered here has a trian- 
gular geometry, the corresponding dual graph possesses a hon- 
eycomb geometry. Note, that we introduced AL extra nodes on 
top and at the bottom of the dual in order to account for the 
free BCs along that direction and to maintain the honeycomb 
structure of the respective graph. This means, a triangular spin 
lattice of size Lx Lis transformed to a honeycomb lattice with 
an over all number of (2L) x (L + 1) nodes. Hence, the topolog- 
ical dual graph associated to the triangular spin lattice is modi- 
fied to some extend. Further note that adjacent extra nodes are 
connected by edges e that carry a weight a>(e) = 0. All other 
edges e on the dual graph get an edge-weight w(e) = JijO-Rjo-Rj. 
Therein, e is assumed to cross a bond Jjj on the spin lattice, 
where Jij couples the two spins ctrj/j, see Fig. |2a). Conse- 
quently, the edge weight on the weighted dual is positive (neg- 
ative), if the corresponding bond on the spin lattice is satisfied 
(broken) with respect to ctr. A pivotal observation is, that there 
exists an equivalence between clusters of adjacent spins on the 
spin lattice that might be flipped in order to decrease the config- 
urational energy of ctr and negative-weighted loops (i.e. closed 
paths) on the weighted dual graph. In this regard, if a loop 
with negative weight on the dual is found, the cluster of spins 
surrounded by this loop can be flipped so as to decrease the 
configurational energy of ctr. Finally, to obtain a GS spin con- 
figuration one needs to find a minimum- weight set of negative- 
weighted loops on the dual graph (see discussion below). This 
set of loops comprises the transition graph for the given real- 
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Figure 2: Illustration of the computation of the transition graph that allows to determine a GS spin configuration on a planar triangular RBIM. The figure illustrates 
a sample system of side length L = 3 and periodic BCs in the horizontal direction, (a) Mapping of the original lattice G (grey edges, triangular geometry), where 
a solid (dashed) line indicates a ferromagnetic (antiferromagnetic) bond, to the weighted dual graph Go (black edges, honeycomb geometry). Note that additional 
edges (dotted lines) where introduced to account for the open boundary conditions in the respective direction. These edges carry zero weight, (b) Minimum weight 
set of loops (bold black edges) on Go as obtained from a mapping to the NWP problem (not shown here, seeRef. 113]). (c) Loop on the dual surrounding a cluster of 
spins on the spin lattice. If the orientation of the spins is chosen as explained in the text, this procedure yields a GS spin configuration. In the figure, spin orientations 
are distinguished by gray filled and non-filled circles. 



ization of the disorder, as illustrated in Fig.|2b). Since initially 
a ferromagnetic reference configuration was chosen, the GS is 
obtained if the orientation of the spins on the spin lattice is cho- 
sen such, that (i) spins within a cluster are aligned in the same 
direction and (ii) spins in adjacent clusters are aligned in oppo- 
site directions. The resulting GS is indicated in Fig. |2c), see 
also Fig. [U 

Here, for the 2d RBIM on a planar triangular lattice, where 
the dual has a honeycomb geometry, the minimum weight set 
of loops on the weighted dual can be obtained by means of a 
mapping to the NWP problem, as explained in [ 1 3.1 . In brief, 
the NWP statement consists in the task to find a minimum- 
weight set of nonintersecting negative-weighted loops for a 
given weighted graph. Therefore, it considers a minimum- 
weight perfect matching problem on an associated auxiliary 
graph with 0{N) edges (provided that the input graph has 0(N) 
edges, as it is the case here), from which the set of loops can 
be deduced. Note that the mapping to the NWP problem yields 
the correct transition graph only for this particular lattice setup, 
since any two-coloring of the spin lattice (i.e. assignment of 
up/down spin orientations) can be composed by loops on the 
dual that do not intersect. That means, each site on the dual 
is an end-node of either or 2 loop segments, as e.g. in Fig. 
|2lb). In contrast, two-colorings of the spins on a square lat- 
tice might involve loops on the dual that involve figure-8 twists. 
That means, each site on the dual is end-node of either 0, 2 or 4 
loop segments. For the latter problem, a different mapping (sl] 
was used recently to obtain exact GSs for 2d ISGs on a square 
lattice with free BCs in at least one direction within polynomial 
time. This mapping was further used to compute "extended" 
GSs for the 2d ISG with fully periodic BCs 1110111 . 

Now, the interpretation of the T - FM-SG transition in 
terms of the NWP problem reads as follows; For small values 
of p, there are only few bonds on the spin lattice that are not 
satisfied by the reference spin configuration. Accordingly, there 
are only few small loops that comprise the transition graph. For 
all nonzero values of p, a sufficiently large lattice will feature 
at least some small loops that surround an elementary plaquette 



on the dual. These small loops correspond to local "manipu- 
lations" of the order parameter (i.e. the magnetization), only. 
Hence, in the thermodynamic limit, the GS has still ferromag- 
netic order (see Fig.flla)). However, if the value of p increases 
and exceeds a critical value pc, large loops appear that have a 
linear extension of the order of the system size and eventually 
span the system along the direction with the periodic boundary 
conditions. These loops represent global manipulations of the 
order parameter, that, in the thermodynamic limit, destroy the 
ferromagnetic order of the GS (see Figs.[Tlb),(c)). 

Once we obtained a GS spin configuration in this manner, we 
compute a MEDW by means of a similar mapping, thoroughly 
explained in Ref. 1 19]. In the following we will use the proce- 
dure outlined above to obtain GSs and to investigate MEDWs 
for the RBIM introduced above. 

3. Results 

As pointed out above, at small values of p there exists an 
ordered ferromagnetic phase, while for large values of p a spin- 
glass ordered phase appears. A proper order parameter to char- 
acterize the respective FM-SG transition is the magnetization 
per spin mi - \ 2, cr,|/L^ for a system of linear extend L. Be- 
low, we perform a finite-size scaling analysis (FSS) in order to 
locate the critical point pc and to estimate the critical exponents 
that describe the scaling behavior of the magnetization in the 



Table 1 : Critical exponents for the 2d RBIM. From left to right: Problem setup 
(SQ=square lattice, TR=fi'iangular lattice), critical exponent of the correlation 
length V, order parameter exponent /?, and exponents 0i and i//[ that characterize 
the scaling of the MEDW excitation energy. The figures for SQ-a are taken from 
Ref. L15J. The figures for SQ-b are taken from Ref. L14J. 
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Figure 3: Results of the finite-size scaling analysis for the binder parameter 
blip), considering dift'erent system sizes L. The main plot shows the unsealed 
data close to the critical point and the insets illustrate the data collapse obtained 
after rescaling the raw data using the scahng assumption discussed in the text 
and scaling parameters as listed in Tab.^ 



vicinity of the critical point. Therefore, we first consider the 
Binder parameter OOll 
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associated with the magnetization. It is expected to scale 
as blip) ~ fi[(p - pc)L^^^], wherein /i[-] signifies a size- 
independent scaling function and v denotes the critical expo- 
nent that describes the divergence of the correlation length as 
the critical point is approached. Here, we simulated triangular 
systems of side length L - 24 through 128 at various values 
of the disorder parameter p. Observables are averaged over 
64000 samples for the largest systems and we used the data 
collapse generated by the scaling assumption above to obtain 
Pc = 0.1584(3)_and v = 1.47(6) with a quahty S = 0.94 of the 
data collapse ll 2lll22ll . see Fig. [3] In general, the above scaling 
relation holds best near the critical point and one can expect that 
there are corrections to scaling off" criticality. As a remedy, we 
restricted the latter scaling analysis to the interval [-0.3, H-0.3], 
enclosing the critical point on the rescaled abscissa. 

Further, the order parameter of the transition is expected to 
scale according to the scaling relation {niiXp)) ~ L'^^^'fjKp - 
Pc)L^^'^], where /2[-] denotes a size-independent function, and 
where the order parameter exponent ^S can be obtained after fix- 
ing V and Pc to the values obtained from the analysis of the 
Binder parameter. The best data collapse (S - 1.01) was ob- 
tained for the choice /? = 0.086(5) (not shown). 

Moreover, an analysis of the average MEDW excitation en- 
ergy (^Zi) according to the scaling assumption ((5£')~L''''/3 [(p- 
Pc)L''"], see Ref. ^, yields the critical point p^ = 0.1586(2), 
in agreement with the above estimate obtained using the Binder 



Figure 4: Results of the finite-size scaling analysis for the average MEDW 
energy (6E). The main plot shows the data collapse after rescaling the raw data 
using the scaling parameters listed in Tab.^and the inset illustrates the unsealed 
data close to the critical point. 



parameter The critical exponents tfri and 0i are listed in Tab. 
[T] Therein, we restricted the scaling analysis to the interval 
[-0.1, -hO.I] and obtained a best data collapse with 5 = 1.63. 

Right at the critical point pc we performed additional simu- 
lations for spin lattices of up to 384 x 384 spins (and 3.6 x 10^ 
samples), i.e. weighted dual graphs of up to 768 x 385 nodes. 
Upon analysis of the data we obtain the estimate /3 - 0.097(8) 
from the scaling behavior of the magnetization, see Fig. |3a). 
We allowed for small deviations from a pure power-law scal- 
ing using a scaling assumption of the form (m) ~ (L + AL)"'^^", 
wherein AL = 0(1). Considering the scaling of the average 
MEDW excitation energy {6E} and using a similar scaling as- 
sumption as above, we found tfri = 0.15(1), see Fig.|5]^b). Both 
these exponents agree within error bars with those obtained ear- 
lier, see Tab.[T] As pointed out above, for the ISG with bimodal 
disorder, there a numerous MEDWs that differ regarding their 
geometric properties. However, here we also analyze the av- 
erage length {{} of the particular MEDWs obtained within the 
simulations, see Fig. |5|b). Therefore, we considered a scal- 
ing according to the form {{} ~ (L + ALY^f, wherein df signi- 
fies the fractal dimension of the MEDWs at pc- We obtained 
df - 1.222(1) (and AL - (9(1)), which is in agreement with 
the value df = 1.222(1) found eariier for the T = FM-SG 



transition for the RBIM on a 2d square lattice, see Ref. 01511 



4. Conclusions 

In the presented article we have illustrated how GSs for the 
2d RBIM on planar triangular lattice graphs can be computed 
by a mapping to the NWP problem. I.e., the problem of find- 
ing a GS spin configuration for a planar 2d triangular RBIM is 
equivalent to the NWP problem on a properly weighted corre- 
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Figure 5: Results of the finite-size scaling analysis at the critical point, where 
the T = FM-SG transition occurs, (a) Scaling behavior of the magnetization 
(m), and, (b) scaHng of the average MEDW length {[) and the MEDW excita- 
tion energy {6E). 



spending dual graph that exhibits a honeycomb structure. Us- 
ing this approach, we have investigated GSs and MEDW exci- 
tations for the respective lattice structure. Therein, a disorder 
parameter could be used to distinguish a ferromagnetic and a 
spin-glass ordered phase. We characterized the corresponding 
r = FM-SG transition by means of a FSS analysis of the mag- 
netization and the MEDW excitation energy. 

In this regard, we found that the values of the critical expo- 
nents obtained here agree within errorbars with those obtained 
earlier for the 2d RBIM on a planar square lattice by consider- 
ing a Gaussian bond distribution with ferromagnetic bias JlSf 
or a bimodal bond distribution [14.1 . as listed in Tab.[T] 

Hence, the results for the triangular lattice structure obtained 
here highlights the universality of the T -Q FM-SG transition. 
Further, note that pc and v found here agree well with the values 
Pc - 0.1583(6) and v - 1.47(9) that characterize the negative- 
weight percolation of loops on 2d lattice graphs with a hon- 
eycomb geometry and fully periodic boundary conditions il3 l. 
Finally, the location of the critical point obtained here via FSS 
analysis is close to the theoretical prediction pcjr - 0.15, that 
was obtained for systems with fully periodic boundary condi- 
tions using the adjoined problem approach i24ll6|] . 
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